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Abstract —  The  use  of  cubic  splines,  instead  of  polynomials,  in  repre¬ 
senting  static  nonlinearities  in  block  structured  models  is  considered.  A 
system  identification  algorithm  for  the  Hammerstein  structure,  a  static- 
nonlinearity  followed  by  a  linear  filter,  is  developed  in  which  the  static 
nonlinearity  is  represented  by  a  cubic  spline.  The  identification  algorithm, 
based  on  a  separable  least  squares  Levenberg-Marquardt  optimization,  is 
used  to  identify  a  Hammerstein  model  of  the  stretch  reflex  EMG  record¬ 
ed  from  a  spinal  cord  injured  patient.  The  resulting  model  provides  more 
accurate  predictions  of  the  reflex  EMG,  even  in  novel  data,  than  more  con¬ 
ventional  models  which  use  polynomial  representations  of  the  nonlineari¬ 
ty.  Furthermore,  the  spline  based  optimization  appears  to  be  less  sensitive 
to  its  initialization  than  a  similar  polynomial-based  approach. 

Keywords —  nonlinear  system  identification,  physiological  model¬ 
ing,  separable  least  squares,  mean  squared  optimization,  Levenberg- 
Marquardt  iteration 


I.  Introduction 

The  Hammerstein  cascade  [10],  a  static  nonlinearity  fol¬ 
lowed  by  a  linear  filter,  is  often  used  to  represent  certain  highly 
nonlinear  systems.  It  is  particularly  useful  for  modeling  sys¬ 
tems,  such  as  the  stretch  reflex,  which  contain  “hard”  nonlin¬ 
earities,  e.g.  rectification  and/or  thresholding.  Indeed,  Ham¬ 
merstein  cascades  can  model  some  nonlinear  systems,  specifi¬ 
cally  those  whose  underlying  structure  is  appropriate,  using  far 
fewer  parameters  than  an  equally  accurate  Volterra  or  Weiner 
series  model. 

Block  structured  models  often  use  polynomials  as  static 
nonlinearities.  Polynomials  are  linear  in  their  parameters,  and 
hence  easy  to  estimate.  Furthermore,  fairly  extreme  nonlinear¬ 
ities  can  be  represented  using  relatively  few  parameters.  Final¬ 
ly,  there  is  a  direct  mathematical  relationship  between  the  poly¬ 
nomial  coefficients  and  corresponding  Volterra  kernels  [10]. 

Several  problems,  however,  have  been  identified  with  poly¬ 
nomial  estimation,  especially  in  regions  where  data  is  sparse 
[3],  Small  disturbances  in  the  data  can  produce  significant 
differences  in  the  interpolated  values.  High  order  polynomial 
solutions  can  also  produce  undesirable  fluctuations,  as  evident 
in  the  polynomial  shown  in  the  upper  panel  of  Fig.  3.  Finally, 
polynomial  extrapolation  is  notoriously  difficult. 

For  these  reasons,  splines  are  often  used  instead  of  poly¬ 
nomials  for  function  approximation  [3],  Indeed,  several  re¬ 
searchers  have  suggested,  but  not  demonstrated,  the  use  of 
splines  in  nonlinear  system  identification  [1],  [2].  Unlike  poly¬ 
nomials,  fitting  splines  requires  a  nonlinear  optimization.  This 
added  difficulty  is  most  likely  the  reason  for  their  limited  appli¬ 
cation  to  date.  Interest,  however,  appears  to  be  growing  in  neu¬ 
ral  networks  applications,  where  the  replacement  of  sigmoidal 
nonlinearities  with  adaptive  splines  is  being  considered  [4]. 

In  this  paper,  we  will  develop  an  identification  algorithm  for 
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Fig.  1.  Block  diagram  of  a  Hammerstein  structure,  consisting  of  a  static  non¬ 
linearity,  followed  by  a  dynamic  linear  system. 


Hammerstein  systems  in  which  the  nonlinearity  is  represent¬ 
ed  by  a  cubic  spline.  The  resulting  technique  will  be  used  to 
identify  models  of  the  stretch  reflex  EMG  from  experimental 
data.  Finally,  identified  models  incorporating  both  polynomial 
and  cubic  spline  nonlinearities  will  be  compared. 

II.  Theory 

The  Hammerstein  model,  shown  in  Fig.  1,  consists  of  a  static 
nonlinearity  followed  by  a  dynamic  linear  system.  Here,  u(t) 
and  y(t)  are  its  input  and  output,  and  x(t)  is  the  signal  between 
the  nonlinear  and  linear  elements. 

The  linear  filter  will  be  represented  by  its  impulse  response 
(IRF),  h(r),  which  will  have  a  memory  length  of  T  samples. 
Its  output  is  computed  using  the  convolution  sum, 

T-r 

V(t)  =  h(r)x{t  -  t)  (1) 

T  — 0 

The  static  nonlinearity,  ra(-),  maps  the  input,  u(t),  to  the 
intermediate  signal,  x(t).  Traditionally,  it  is  modeled  with  an 
order  Q  polynomial,  in  which  case, 

Q 

x(t )  =  m{u{t))  =  Y,cqu\t)  (2) 

g=o 

where  cq  is  the  </’th  order  polynomial  coefficient.  Regardless 
of  how  the  nonlinearity  is  represented,  the  output  of  a  Ham¬ 
merstein  cascade  can  be  written  as, 

T-l 

y(t )  =  ^2  -  T))  (3) 

T— 0 

In  the  special  case  of  a  polynomial  Hammerstein  model,  where 
m(-)  is  parameterized  with  a  polynomial,  (3)  becomes, 

T-l  Q 

y (t)  =  cqh(T)u9(t  -  r)  (4) 

t— 0  q— 0 


Report  Documentation  Page 

Report  Date 

Report  Type 

Dates  Covered  (from..,  to) 

250CT2001 

N/A 

- 

Title  and  Subtitle 

Contract  Number 

Identification  of  a  Hammerstein  Model  of  the  Stretch  Reflex 

EMG  using  Cubic  Splines 

Grant  Number 

Program  Element  Number 

Author(s) 

Project  Number 

Task  Number 


Work  Unit  Number 

Performing  Organization  Name(s)  and  Address(es)  Performing  Organization  Report  Number 

Dept.  Elec.  &  Comp.  Eng.,  Univ.  Calgary,  2500  University  Dr. 

NW,  Calgary,  AB,  T2N  1N4,  Canada" 

Sponsoring/Monitoring  Agency  Name(s)  and  Address(es)  Sponsor/Monitor’s  Acronym(s) 

US  Army  Research  Development  &  Standardization  Group  - 

(UK)  PSC  803  Box  15  FPO  AE  09499-1500  Sponsor/Monitor’s  Report  Number(s) 

Distribution/ Availability  Statement 

Approved  for  public  release,  distribution  unlimited 

Supplementary  Notes 

Papers  from  the  23rd  Annual  International  conference  of  the  IEEE  Engineering  in  Medicine  and  Biology  Society, 
October  25-28,  2001,  held  in  Istanbul,  Turkey.  See  also  ADM001351  for  entire  conference  on  cd-rom. 

Abstract 


Subject  Terms 


Number  of  Pages 

4 


2 


Techniques  based  on  separable  least  squares  (SLS)  opti¬ 
mizations  [9]  have  been  developed  for  the  identification  of 
polynomial  Hammerstein  models  [11],  Briefly,  the  Hammer- 
stein  cascade  is  represented  by  a  parameter  vector  containing 
the  filter  weights  and  polynomial  coefficients, 

9  =  [  h{ 0)  . . .  h(T  -  1)  I  Co  ...  CQ  ]  T  (5) 

=  [oJ\el]T 

Its  output  is  written  y(t,  0).  explicitly  showing  the  dependence 
on  the  parameter  vector.  Thus,  y(t,  9)  is  computed  using  (4), 
where  the  IRF  and  polynomial  coefficients  are  taken  from  9, 
defined  in  (5).  Furthermore,  identifying  the  model  is  equiva¬ 
lent  to  finding  the  parameter  vector  that  minimizes  the  MSE, 

1  N 

Vn{0)  =  (S/W  “  ?/(*> 9 ))2 


The  SLS  identification  algorithm  [11]  hinges  on  the  obser¬ 
vation  that  for  any  given  choice  of  polynomial  coefficients,  the 
output  (4)  is  a  linear  function  of  the  filter  weights.  This  is 
reflected  in  (5),  where  9  is  divided  into  linear  and  nonlinear 
parameters,  9i  and  9n  respectively.  Since  (4)  is  linear  in  9i,  its 
optimal  value,  corresponding  to  any  choice  of  polynomial  co¬ 
efficients,  9n,  can  be  found  in  closed  form  by  solving  a  linear 
regression.  Thus,  9i  is  a  function  of  9n,  written  9i(9n).  Sim¬ 
ilarly  y(9n ),  and  Vj v(6n)  are  also  functions  of  the  nonlinear 
parameters  alone.  Thus,  an  iterative,  nonlinear,  optimization  is 
only  needed  to  find  9n.  In  this  paper,  we  use  the  Levenberg- 
Marquardt  algorithm.  Thus,  9n  is  updated  using: 

€+1)  =  +  (Jj  Js  +  Skl)~ 1  Jj  e  (6) 

where  e  is  a  vector  whose  f’th  element  is  the  error  e(t)  = 
y(t)  —  y(t),  I  is  an  identity  matrix,  and  6  k  is  a  regularization 
parameter  used  to  control  the  convergence  rate  and  stability. 
Js ,  the  Jacobian,  is  matrix  whose  [f ,  to]  entry  contains  the  par¬ 
tial  derivative  of  the  model  output  at  time  t  with  respect  to  the 
?n’th  element  of  9n.  Thus, 


Js(t,  to) 


dy(t) 

dOn(m) 


(7) 


Note,  however,  that  in  computing  Js.  the  dependence  of  0/,  the 
linear  parameters,  on  9n  must  be  taken  into  consideration.  This 
can  be  done  as  follows.  First  compute  <//  and  ./„,  the  Jacobians 
with  respect  to  the  linear  and  nonlinear  parameters,  assuming 
that  the  parameters  are  independent  of  each  other.  Then,  it  can 
be  shown  [9]  that 

Js  =  (I-  P,)Jn  (8) 

where  Pi  =  Ji(J^  Ji)~x  Jf  is  an  orthogonal  projection  onto 
the  columns  of  the  linear  Jacobian, 

Thus,  to  use  the  SLS  algorithm,  compute  the  Jacobian  in  the 
usual  way,  and  then  separate  it  into  linear  and  nonlinear  parts. 
Insert  these  into  (8)  to  compute  Js,  used  in  the  L-M  step  (6). 


A.  Use  of  Cubic  Splines 

The  SLS  algorithm  depends  on  the  output,  (4),  being  a  linear 
function  of  the  IRF  weights,  but  does  not  depend  on  the  rep¬ 
resentation  of  the  static  nonlinearity,  to(-).  All  that  is  needed 
to  use  a  different  nonlinearity,  (3),  is  to  compute  the  nonlinear 
Jacobian,  with  respect  to  its  parameters. 

In  this  paper,  the  static  nonlinearity,  x  =  m(u),  will  be 
modeled  by  a  cubic  spline  instead  of  a  polynomial.  A  cu¬ 
bic  spline  is  defined  by  a  series  of  M  knots ,  ( Uj,Xj ),  for 
j  =  1,2, ...  ,M,  with  u\  <  U2  <  ...  <  Um-  Thus,  the 
parameter  vector  describing  the  spline  can  be  formed, 

9n  =  [  Ml  ...  Um  Xl  ...  xM  } 1  (9) 

In  between  each  pair  of  knots,  the  spline  is  defined  by  a  third 
degree  polynomial.  These  are  chosen  so  that  the  spline  and 
its  first  two  derivatives  are  all  continuous  functions.  Hence, 
let  u  be  a  point  in  the  input,  between  knots  j  and  j  +  1,  (i.e. 
Uj  <  u  <  Uj+ 1).  The  corresponding  point  in  the  output  of  the 
spline  is  computed  from  [8], 

x  =  A(u)xj  +  B{u)Xj+ 1  +  C(u)Xj  +  D(u)x'j+1  (10) 

where  A(u),  B(u),  C(u)  and  D(u)  are  defined  by 


A(u)  =  U°+l  U 
uj+ i  uj 

B(u )  =  1  —  A(u) 

C{u)  =  ^(A*{u)~  A{u)){uj+1-ujf 

D(u )  =  ^(H3(w)  -  B(u))(uj+ 1  -  uj)2 

Note,  however,  that  (10)  depends  on  the  second  derivative  of 
the  the  spline,  x"  =  d2x/du2,  at  the  two  adjacent  knots.  At 
first  glance,  it  would  appear  that  these  second  derivatives  are 
required  to  define  the  spline  completely.  However,  by  forcing 
the  first  derivative  of  the  spline  to  be  continuous  across  the 
knots,  one  arrives  at  the  equation  [8], 


ui~ui- 1 

6  Xj~  1 


j  +  i-ttj-i  m// 

3  xj 


i+l~x 
Uj  +  1~  Uj 


+i  ~uj 

6  xj+ 1 

X-i  —Xj 


(11) 


for  j  =  2, . .  .,M  —  1.  Thus,  we  have  M  —  2  equations 
with  which  to  define  x"  at  M  locations,  which  is  an  under¬ 
determined  problem.  Typically,  this  is  resolved  by  arbitrari¬ 
ly  choosing  the  second  derivatives  at  the  two  end  points.  In 
this  paper,  natural  splines  will  be  employed,  where  the  second 
derivative  is  set  to  zero  at  the  first  and  last  knots.  Thus  (1 1)  can 
be  used  to  solve  for  x"  at  the  remaining  M  —  2  knots.  Once 
these  have  been  computed,  (10)  can  be  used  to  compute  the 
output  of  the  cubic  spline. 

To  use  a  cubic  spline,  instead  of  a  polynomial,  in  the  SLS 
Hammerstein  identification  algorithm,  one  need  only  compute 
the  appropriate  nonlinear  Jacobian, 


Jn 


dx(t) 
du\  ’ 


dx(t)  dx(t ) 
8um  ’  dxi  ’ 


dx(t ) 

dxM 


(12) 


and  then  insert  it  into  (8)  to  compute  the  separated  Jacobian  in 
the  Levenberg-Marquardt  step  (6).  The  rest  of  the  algorithm  is 
used  without  modification. 
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Typical  Experimental  Data 


Time  (s) 

Fig.  2.  Extract  from  a  typical  experimental  trial  showing  4  seconds  of  position, 
velocity  and  GS  EMG. 

III.  Experimental  Results 

To  evaluate  the  utility  of  the  cubic  spline  algorithm,  we  used 
it  to  estimate  stretch  reflex  dynamics,  the  relationship  between 
the  ankle  velocity  and  the  resulting  EMG,  in  spastic,  spinal 
cord  injured  (SCI)  patients.  Previously,  this  system  has  been 
modeled  as  a  Hammerstein  cascade  [5]  in  which  the  static  non¬ 
linear  element  was  found  to  resemble  a  half-wave  rectifier  and 
the  IRF  of  the  linear  element  resembled  a  delayed  impulse. 
Since  this  model  is  being  investigated  for  potential  clinical  use 
[7],  methods  that  find  efficient,  unbiased  estimates  of  its  ele¬ 
ments  may  prove  to  be  very  significant. 

The  experimental  procedures  have  been  described  in  detail 
elsewhere  [6],  In  summary,  the  subjects,  due  to  a  previous 
spinal  cord  injury,  demonstrated  an  incomplete  loss  of  motor 
function,  spasticity  and  hyper-active  stretch  reflexes.  They  lay 
on  their  backs  with  their  left  foot  attached  to  a  rotary  hydraulic 
actuator  by  a  custom  fitted  fiber-glass  boot.  Ankle  torque  and 
position,  and  the  EMG  over  the  Gastrocnemius-Soleus(GS) 
muscles  were  recorded.  A  broad-band  position  perturbation, 
whose  spectrum  was  shaped  to  preserve  the  stretch  reflex  [6], 
was  applied  while  the  subject  maintained  a  constant  back¬ 
ground  contraction.  Figure  2  shows  four  of  the  30  seconds, 
sampled  at  200  Hz,  of  position,  velocity  and  EMG  data  used  in 
the  analysis.  Note  that  the  velocity  was  calculated  by  numeri¬ 
cally  differentiating  the  measured  position. 

The  technique  developed  in  [11]  was  used  to  fit  a  polyno¬ 
mial  Hammerstein  model  between  the  first  5000  points  of  the 
velocity  and  EMG  records,  hereafter  called  the  training  data. 
The  starting  point  for  this  SLS  optimization  was  the  nonlinear¬ 
ity  identified  using  a  conventional  correlation  based  scheme 
[5].  The  resulting  model,  shown  in  Fig.  3  with  dash-dotted 
lines,  was  then  used  to  predict  the  remaining  1000  points  of 
EMG  data,  the  validation  segment.  The  polynomial  Hammer¬ 
stein  model  accounted  for  95.2%  of  the  signal  variance  in  the 
training  sample  and  94.5%  in  the  validation  set. 

Next,  a  cubic  spline  Hammerstein  cascade  was  identified 


Static  Nonlinearity 
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Fig.  3.  Hammerstein  models  of  the  stretch  reflex  EMG  using  polynomial  and 

cubic  spline  nonlinearities  identified  from  experimental  data. 

from  the  training  data.  The  initial  spline  had  5  knots  along 
the  line  m(u)  =  u,  equally  spaced  over  the  range  of  the  exper¬ 
imental  input.  The  SLS  iteration  described  in  this  paper  was 
then  used  to  search  for  the  optimal  cubic  spline  nonlinearity 
and  corresponding  linear  IRF.  The  results  of  this  identification, 
referred  to  as  “spline  one”,  are  shown  as  solid  lines  in  Fig.  3. 
The  nonlinearity  is  shown  in  the  upper  panel:  the  knots  are 
shown  as  circles,  with  the  behavior  between  knots  shown  as  a 
solid  line.  The  prediction  accuracy  of  this  model  was  95.3% 
VAF  in  the  training  data,  and  95.6%  VAF  in  the  validation  sam¬ 
ple.  Therefore,  using  the  cubic  spline  algorithm  decreased  the 
residual  variance  by  2.1%  and  25.0%  in  the  training  and  vali¬ 
dation  sets,  respectively. 

If  a  general  model  of  the  system  has  been  previously  devel¬ 
oped,  points  for  the  initial  spline  can  be  taken  from  its  nonlin¬ 
earity.  For  example,  a  cubic  spline  model,  referred  to  as  “spline 
two”,  was  generated  from  the  same  initial  model  as  the  poly¬ 
nomial  cascade.  The  results  are  shown  in  Fig.  3  as  dashed  lines 
with  the  spline  knots  represented  by  asterisks.  Spline  two  ac¬ 
counted  for  95.5%  of  the  variance  in  the  training  set  and  96.0% 
in  the  validation  set,  giving  a  decrease  of  6.7%  and  37.5%  of 
the  residual  variance  in  the  two  data  segments  with  respect  to 
the  polynomial-based  cascade. 

Comparing  the  models  produced  by  the  two  algorithms  il¬ 
lustrates  several  potential  advantages  of  the  spline-based  ap¬ 
proach.  First,  its  models  produced  more  accurate  predictions 
than  did  the  polynomial-based  scheme.  Secondly,  the  “spline 
one”  model  resulted  from  an  initialization,  a  simple  linear  fit, 
that  contained  no  a  priori  information  about  the  shape  of  the 
nonlinearity,  suggesting  that  the  spline-based  algorithm  may 
be  less  sensitive  to  its  initialization  than  the  polynomial  ap¬ 
proach.  Furthermore,  neither  spline -based  model  included  the 
oscillations  at  negative  velocities  present  in  the  polynomial 
curve  (see  Fig.  3  upper  panel).  Finally,  the  dynamic  linear 
subsystem  (lower  panel)  of  the  “spline  two”  model  shows  a 
steeper  response  with  a  longer  rest  time  at  zero,  which  is  more 
consistent  with  the  propagation  delay  expected  in  the  reflex. 
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Fig.  4.  Upper  panel:  input-output  scatter  plot  of  spline  two.  Lower  panel: 

probability  density  histogram  of  input  velocity. 

In  both  the  polynomial  and  spline  nonlinearities,  there  is  a 
large  negative  deflection  between  input  velocities  of  0.6  and 
1.1  rad/sec.  However,  this  deflection  appears  to  be  well  sup¬ 
ported  by  the  experimental  data.  Indeed,  the  upper  panel  of 
Fig.  4  shows  a  scatter  plot  of  the  input  and  output  of  the  non¬ 
linearity  from  spline  two.  Since  much  of  the  deflection  ap¬ 
pears  as  a  series  of  solid  lines,  it  is  evidently  well  supported  by 
the  training  data.  The  lower  panel  shows  a  probability  density 
histogram,  also  estimated  from  the  training  data,  for  the  non¬ 
linearity  input.  This  suggests  that,  except  near  zero  velocity, 
all  points  in  the  input  are  probed  with  almost  equal  frequen¬ 
cy.  Since  the  parameters  describing  a  cubic  spline  (i.e.  the 
knot  positions)  have  primarily  local  effects  on  its  shape,  it  is 
unlikely  that  the  deflection  is  an  artifact  due  to  the  analysis. 

IV.  Discussion 

In  this  paper,  we  have  developed  a  cubic  spline  identifica¬ 
tion  method  for  Hammerstein  systems  and  used  it  to  identify 
models  of  the  stretch  reflex  EMG  in  SCI  patients.  Both  spline- 
based  models  developed  more  accurate  predictions  than  a  poly¬ 
nomial  based  cascade  identified  from  the  same  data.  This  ap¬ 
plied  to  both  the  training  and  validation  samples. 

Iterative  optimizations,  such  as  the  cubic-spline  based  tech¬ 
nique  developed  in  this  study,  are  prone  to  finding  sub-optimal 
local  minima.  For  example,  the  results  obtained  by  the  poly¬ 
nomial  method  [11]  are  dependent  on  the  choice  of  initial 
polynomial  coefficients.  In  this  study,  see  Fig. 3,  changing 
the  initial  spline  knots  produced  slight  variations  in  the  final 
Hammerstein  cascade.  During  the  data  analysis,  several  other 
initial  splines  were  tested  (results  not  shown).  Starting  with 
m(u)  =  k,  a  constant,  led  to  a  model  in  which  the  nonlinearity 
was  offset  and  the  IRF  was  highly  oscillatory.  This  model  was 
less  accurate  than  even  the  polynomial  cascade.  On  the  other 
hand,  when  the  initial  knots  were  taken  from  m(u)  =  u2  or 
m(u)  =  u3,  the  optimization  found  models  similar  to  splines 
one  and  two.  Overall,  the  cubic  spline  method  appeared  to  be 


less  dependent  on  its  initial  parameters  than  the  polynomial- 
based  algorithm.  Future  research  will  analyze  the  convergence 
behavior  of  the  cubic  spline  Hammerstein  algorithm. 

The  price  of  the  cubic  spline -based  algorithm  is  an  increased 
number  of  parameters  needed  to  characterize  the  nonlinearity. 
A  cubic  spline  Hammerstein  cascade  is  still  much  more  effi¬ 
cient  than  an  equivalent  Volterra  series  model,  but  uses  more 
parameters  than  a  polynomial  Hammerstein  model.  In  this 
study,  a  cubic  spline  with  five  knots,  which  has  ten  param¬ 
eters,  was  compared  to  a  sixth  order  polynomial,  with  seven 
coefficients.  With  the  cubic  spline  initializations,  including  ad¬ 
ditional  knots  provided  only  a  small  increase  in  accuracy,  but 
beginning  with  fewer  than  five  resulted  in  significantly  inferior 
models.  Similarly,  increasing  the  nonlinearity  order  had  little 
effect  on  the  accuracy  of  polynomial-based  models.  Finally, 
increasing  the  memory  length  of  the  IRFs  led  to  slightly  small¬ 
er  residuals  for  both  the  cubic  spline  and  the  polynomial-based 
cascades. 

These  results  suggest  that  cubic  splines  may  be  used  in  place 
of  polynomials  in  a  wide  variety  of  block-structured  models 
currently  used  to  represent  physiological  systems.  It  appears 
that  iterative  identification  algorithms  based  on  splines  may 
be  less  dependent  on  their  initialization  than  similar  poly¬ 
nomial-based  schemes.  Furthermore,  models  based  on  splines 
are  likely  to  be  relatively  immune  to  the  severe  interpola¬ 
tion/extrapolation  problems  associated  with  polynomials. 
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